/*!  \file  SteVR.cpp
 *	  \brief  
 * 
 *   \details   details
 *   \author    Yordan Madzhunkov German Research School for Simulation Sciences GmbH
 *   \version   1.0
 *   \date      07/07/2013 11:17:10 AM
 *   \copyright GNU Public License.
 */
#include "SteVR.h"
#include "Matrix.h"
#include "setThreads.h"
extern "C" {
void dstevr_( 
const char   *  jobz, 
const char   * range, 
const int    *     n, 
double *     d, 
double *     e, 
const double * vl, 
const double * vu, 
const int    * il, 
const int    * iu, 
const double * abstol,
int    * m, 
double * w, 
double * z, 
const int    * ldz, 
int    * isuppz, 
double * work, 
const int    * lwork, 
int    * iwork,
const int    * liwork, 
int    * info );
}

SteVR::SteVR(): TriEigen("Lapack dstevr") {	
	setCompileTime(__DATE__, __TIME__); 
	canBeParallel = false;
}
bool SteVR::setNumberOfThreads(int nThreads_) {	
	assert(canBeParallel);
	SET_NUM_THREADS(nThreads_);
	nThreads = nThreads;
	return true;
}
#include <iostream>
using namespace std;
void SteVR::compute(const Matrix & T, Matrix & Sigma, Matrix & Q ){
	reset();
	assert(T.isTriDiagonal() && Sigma.isDiagonal() && Q.isOrthogonal());
	char jobz = 'V', range = 'A'; int n = T.cols(), info = -1, ldq = Q.lda();
	double * e = NULL; allocate(e, n); 
	double * d = NULL; allocate(d, n); 
	int * isuppz = NULL; allocate(isuppz, 2 * n);
	memcpy(d, T.getData()    ,    n * sizeof(*d));
	memcpy(e, T.getOffDiagonal(), (n-1)* sizeof(*d));
	double vl = 0; double vu = 1; 
	int il = 1; int iu = n; int m;
	double abstol = 0;  
	double * w = Sigma.getModData();
	double * q = Q.getModData(); Q.makeColMajor();
	int lwork = -1, liwork = -1;
	double wopt; int iwopt;
	dstevr_(&jobz, &range, &n, d, e, &vl, &vu, &il, &iu, &abstol, &m, w, q, 
		&ldq,	isuppz, &wopt, &lwork, &iwopt, &liwork,	&info);
	success = false;
	if (info == 0) {
		lwork = (int) wopt; liwork = iwopt;
		double *  work = NULL; allocate( work, lwork);
		int    * iwork = NULL; allocate(iwork, liwork);
		dstevr_(&jobz, &range, &n, d, e, &vl, &vu, &il, &iu, &abstol, &m, w, q, 
			&ldq, isuppz, work, &lwork, iwork, &liwork,	&info);
		success = info == 0;
		delocate( work);
		delocate(iwork);
	}
	delocate(isuppz);
	delocate(d); 
	delocate(e); 
};
BiNESteVR::BiNESteVR() : BiSVD( "NE dstevr" ) {
	aUpper = true; 
	aLower = false; 
	aUColMajor = true; 
	aVColMajor = true; 
	VColVec = false;
	setTriEigen(*(new SteVR()));
	setCompileTime(__DATE__, __TIME__); 
};
BiNE2SteVR::BiNE2SteVR() : BiSVD( "NE2 dstevr" ) {
	aUpper = true; 
	aLower = false; 
	aUColMajor = true; 
	aVColMajor = true; 
	VColVec = false;
	setTriEigen(*(new SteVR()));
	setCompileTime(__DATE__, __TIME__);
};
//REGISTER(SteVR);
//REGISTER(BiNESteVR);
//REGISTER(BiNE2SteVR);
